#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sun Jun  8 15:19:35 2025

@author: fenghongli
"""

import os
import pandas as pd
import numpy as np
from statsmodels.tsa.api import VAR
import matplotlib.pyplot as plt


input_dir = 'garch_volatility'
vol_dict = {}

# 读取所有volatility文件
for file in os.listdir(input_dir):
    if file.endswith('_volatility.csv'):
        ts_code = file.replace('_volatility.csv', '')
        df = pd.read_csv(os.path.join(input_dir, file))
        df['trade_date'] = pd.to_datetime(df['trade_date'], format='%Y%m%d')
        df.set_index('trade_date', inplace=True)
        vol_dict[ts_code] = df['volatility']

# 合并
vol_df = pd.DataFrame(vol_dict)

# ⚠️ 修复关键点：使用线性插值补全缺失值（保留所有资产）
vol_df = vol_df.interpolate(method='linear', limit_direction='both')

# ⚠️ 再删除极少量仍缺失的（若有）
vol_df = vol_df.dropna()

# 保存
vol_df.to_csv('merged_volatility_panel.csv')
print(f"合并完成，资产个数: {vol_df.shape[1]}，有效日期数: {vol_df.shape[0]}")



# 读取面板数据
df = pd.read_csv('merged_volatility_panel.csv', index_col='trade_date', parse_dates=True)

window_size = 200  # 滚动窗口长度
horizon = 10       # 预测步长
total_so_list = [] # 总溢出指数
dates = []

# 溢出指数计算函数（简化版）
def compute_dy_so_normalized(var_model, steps):
    """
    使用标准化方差分解矩阵计算 Diebold-Yilmaz 总溢出指数
    """
    fevd_matrix = var_model.fevd(steps).decomp[-1]  # T期预测误差方差分解矩阵（N×N）
    
    # 标准化每一行（按行归一化）
    fevd_normalized = fevd_matrix / fevd_matrix.sum(axis=1, keepdims=True)
    
    # 剔除对角线（自身影响）
    np.fill_diagonal(fevd_normalized, 0)
    
    # 总溢出指数 = 所有非对角元素均值 × 100
    N = fevd_normalized.shape[0]
    total_spillover = 100 * fevd_normalized.sum() / N
    
    return total_spillover


# 滚动窗口
for i in range(window_size, len(df)):
    window_data = df.iloc[i - window_size:i]
    try:
        model = VAR(window_data)
        result = model.fit(maxlags=1)
        so = compute_dy_so_normalized(result, steps=horizon)
        total_so_list.append(so)
        dates.append(df.index[i])
    except Exception as e:
        total_so_list.append(np.nan)
        dates.append(df.index[i])
        print(f"窗口{i} VAR计算失败：{e}")

# 画出总溢出指数时序图
plt.figure(figsize=(12, 5))
plt.plot(dates, total_so_list, label='总溢出指数（DY）', color='darkred')
plt.title('DY Total Spillover Index (Rolling VAR)')
plt.xlabel('time')
plt.ylabel('Spillover Intensity(%)')
plt.grid(True)
plt.tight_layout()
plt.savefig('total_dy_spillover_index.png', dpi=300)
plt.show()
